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The  physical  procaa—  of  tha  aarth'a  ataoaphare  can  bo 
sods lied  toy  a  systaa  of  hydrodynaslc  aquations.  This  systaa  of 
aquations  cannot  too  solved  directly  unless  eany  siaplif ying  as- 
eusp t ions  are  aada,  savsraly  Halting  haw  realistically  tha 
actual  ataospherlc  procsssss  are  sieulatod.  In  order  to  produce 
accurate  weather  forecasts  the  full  systse  of  aquations  must  be 
solved  in  four  di pensions.  In  practice  global  weather  forecasts 
are  produced  at  various  eetsorological  canters  around  the  world 
by  treating  these  equations  as  an  initial  value  problsa  and 
integrating  forward  in  ties  to  produce  forecasts.  The  solution 
of  this  problsa  requires  the  use  of  advanced  vector  processors  as 
the  nuaber  of  coaputations  involved  is  staggering.  The  forecast¬ 
ing  protoloe  was  foraulatad  quite  succinctly  by  the  Norwegian 
pioneer  in  weather  forecasting*  V.  Bjerknes  C13,  when  he  defined 
the  necessary  and  sufficient  conditions  for  a  successful  systea 
in  an  article  written  in  1922  to  bei 

i)  A  sufficiently  accurate  knowledge  of  the  state  of  the 
ataosphere  at  the  initial  tlae. 

ii)  A  sufficiently  accurate  knowledge  of  the  laws  according 
to  which  one  state  of  the  ataoaphare  develops  froe  another. 

Bjerknes'  discover les  of  the  hydrodynaaical  nature  of  the 
weather  problsa  led  several  European  nations*  especially  the 
Scandinavian  countries*  to  begin  collecting  observations  of  the 
state  of  the  ataosphere.  This  data  collection  led  L.  F. 

Richardson  C23  to  try  describing  initial  conditions  froa  a  hand 

'  >  . 

analysis  and  projecting  the  state  of  the  ataosphere  to  the  future 
froa  the  hydrodynaaical  equations.  The  task  was  monumental  as 


Richardson  MtiMtad  that  a  warehouse  of  64,000  paopla  using  tha 
— chanieal  calculator  of  tha  day  could  just  foracast  tha  stata  of 
tha  atmosphere  at  tha  rata  that  it  was  actually  evolving.  Unfor¬ 
tunately,  aany  factors,  discovered  later  in  tha  1940's,  kept 
Richardson  froa  asking  a  successful  forecast. 

Tha  aagnituda  of  tha  weather  forecasting  pr obi aa  required 
tha  davalopaant  of  electronic  computers  for  even  siaple 
solutions.  Tha  Electronic  Nueerical  Integrator  and  Coaputer 
(ENIAC)  developed  in  the  late  1940's  allowed  Charney,  Fjortoft, 
and  Von  Neuaann  C31  to  succeed  in  asking  a  reasonable  24  hour 
forecast.  Their  hydrodynaaical  model  was  siaplified  to  filter 
gravity  waves  while  allowing  weather  patterns  to  develop  in  a 
aanner  siailar  to  that  observed  in  the  atmosphere.  Their  initial 
conditions  of  pressure  heights  were  derived  by  hand  and  the 
result  gridded  and  typed  into  the  computer. 

With  the  rapid  development  of  coaputers  over  the  past  thirty 
years,  it  has  become  possible  to  use  numerical  techniques  to 
integrate  the  full  set  of  hydrodynamic  equations  forward  in  time 
to  produce  improved  weather  forecasts.  As  V.  Bjerknes  predicted, 
accurate  forecasts  require  more  than  just  accurate  treatment  of 
the  physical  processes  of  the  atmosphere,  they  also  require 
accurate  specifications  of  the  initial  state  froa  non-uni f oral y 
located  observations.  Panofsky  C43,  Bergthorsson  and  DOfli  [5], 
and  Cressaan  C63  pioneered  methods  to  use  the  computer  to  obtain 
a  weather  analysis  froa  observational  data.  This  process  of 
combining  observation  values  with  a  background  field  is  called 
objective  analysis.  For  the  most  part,  these  original  objective 
analysis  techniques  were  weighted  average  schemes  that  depended 


upon  propor  pMKiliutian  of  Mv«r«l  pariMtari,  usually  obtainad 
la  an  ad  hoe  way.  Today,  aoat  of  tha  world's  woathsr  csntsrs  usa 
statistical  abjective  analysis  tschnlquss  bassd  on  ths  work  of 
•andin  C73  to  provids  Initial  conditions  for  their  atmospheric 
forecast  nodal s. 

In  practice  two  sourcas  of  information  arm  combi nad  to 
produce  an  objective  analyaisi  observations  of  atmospheric 
variables  and  a  forecast  made  by  the  atmospheric  model  from  a 
previous  analysis.  The  forecast  is  commonly  called  the  'first 
guess*  to  the  analysis  or  the  'background*.  Because  the  forecast 
is  hardly  a  guess,  the  term  'forecast  background*  is  used  in  this 
paper  to  emphasize  that  the  background  to  the  analysis  was 
derived  from  a  forecast  made  earlier.  The  observations  of 
temperature,  wind,  and  moisture  are  made  by  in  situ  instruments 
attached  to  balloons,  aircraft,  and  ships  or  from  remote  instru¬ 
ments  aboard  satellites  or  on  the  earth's  surface.  The  result  is 
a  collection  of  observations  of  varying  degrees  of  accuracy  taken 
at  various  times.  The  statistical  analysis  schemes  have  been 
designed  to  'optimally'  combine  these  observations  with  the  fore¬ 
cast  background  to  produce  the  initial  conditions  required  by  the 
numerical  forecast  model.  The  optimality  of  these  schemes 
directly  depends  upon  how  well  the  statistical  properties  of  the 
errors  of  the  forecast  background  and  the  observations  are 
defined.  In  practice,  the  schemes  are  multivariate  in  the  sense 
that  they  are  used  to  simultaneously  analyze  multiple  related 
dependent  variables  from  measured  values. 

In  this  paper  we  will  deal  with  the  representation  of  the 
statistics  of  the  forecast  background  error.  In  particular,  the 


■octal ling  of  the  spatial  autocovar i ance  of  ths  srror  for  ths 
primary  variabls  is  exaainad.  Early  varsions  of  Bandin's  osthod 
usad  a  si apis  exponential  function  to  aodal  ths  autocovariancs. 
Although  this  aodal  is  simple,  it  failad  to  bs  sufficiently 
flexible  to  dascribs  details  of  the  statistics  derived  froa 
actual  data.  A  search  of  the  literature  revealed  that  many 
models  are  available,  but  tests  of  their  abilities  in  fitting 
background  statistics  for  an  actual  forecast  model  have  not  been 
conducted.  The  mathematical  and  precision  limitations  of  various 
models  have  been  determined  and  are  described  in  this  paper. 

Optimum  interpolation  (01),  which  is  sometimes  more  prac¬ 
tically  referred  to  as  statistical  interpolation  (SI),  is  applied 
to  compute  the  corrections  to  the  background  field.  This  is  done 
by  first  interpolating  the  background  to  the  non-uni formly  lo¬ 
cated  observation  locations,  and  then  computing  the  difference 
between  the  observed  and  background  value.  If  observations  were 
exact,  this  would  be  the  background  error  measured  at  discrete 
locations.  These  measurements  of  background  error  are  then  used 
by  01  to  compute  a  correction  field  on  a  uniform  grid,  which  is 
added  to  the  background  to  produce  the  analysis. 

The  full  development  of  the  equations  for  a  multivariate 
application  of  01  is  given  in  several  papers,  including 
Rutherford  C83,  Schlatter  C93,  Schlatter,  et.  al .  C101,  Bergman 
C113,  Lorenc  C12D,  Thldbaux  C13D,  and  Thidbaux  and  Pedder  C143. 

A  brief  outline  of  the  method  is  given  in  the  following.  For  a 
collection  of  estimates  of  the  error  at  scattered  points,  it  is 
desired  to  estimate  the  value  of  the  error  at  the  grid  points. 

01  approximates  these  values  by  a  linear  combination  of  the  known 


values  dallirad  so  that  tha  aw  pact  ad  aaan  squared  arror  ovar  aaaa 
anaaabla  of  realisations  la  minimized.  Thia  raquirea  that  tha 
atatlatical  propartiaa  (covariancaa  between  variables)  ba  known. 
Stationari ty  (lndapandanca  of  tha  particular  grid  point)  of  tha 
atatiatical  paraaatara  ia  raquirad  for  a  tractabla  problaa.  Tha 
weights  uaad  in  tha  linaar  combination  are  obtainad  from  tha 
solution  of  a  car tain  aystam  of  linaar  aquations,  tha  coafficlant 
matrix  baing  tha  matrix  of  covariancaa  batwaan  tha  background 
plus  obaarvad  arrors  at  tha  observation  points.  Tha  positive 
definiteness  of  this  matrix  plays  an  Important  role,  both  the¬ 
oretically  and  computationally. 

A  discussion  of  multivariate  covarianca  functions,  proper — 
ties  they  must  satisfy,  and  methods  of  obtaining  such  functions 
are  discussed  in  section  2.  Experiments  with  savaral  families  of 
covarianca  functions  in  fitting  background  arror  statistics  and 
tha  resulting  performance  in  a  statistical  interpolation  schema 
are  described  in  section  3.  Section  4  summarizes  tha  results 
and  suggests  soma  future  work. 

2.0  iiulti variata  Covarianca  Functions 

2. 1  Oanaral  Development 

Tha  theory  of  covarianca  functions  and  that  of  positive 
definite  functions  go  hand-in-hand.  Positive  definiteness  of 
matrices  such  as  occur  in  our  application  are  equivalent  to  the 
spatial  covariance  function  for  the  background  errors  being  posi¬ 
tive  definite.  Positive  definite  functions  are  characteri zed  by 
Bochner’s  Theorem  C153,  which  states  that  a  function  is  positive 
(semi ) definite  if  and  only  if  its  Fourier  transform  is 


wonwf  tlvt.  Altariutivaly,  the  covariance  function  la  ths 
Fourier  (cosino)  transform  of  a  probability  density  (nonnegative) 
function,  because  of  the  application,  our  interest  is  in  posi¬ 
tive  definite  functions  that  are  s sooth  in  the  senes  that  certain 
partial  derivatives  exist.  An  excellent  reference  for  positive 
definite  functions  is  Stewart  C163. 

For  completeness,  a  derivation  of  the  covariance  functions 
for  variables  related  through  differentiation  is  given  hare. 
Suppose  that  it  is  wished  to  analyse  three  related  depandent 
variables,  requiring  that  the  corrections  obtained  via  01  (or 
■ore  correctly,  81)  will  not  upset  the  relationship  between  the 
predicted  values  of  the  variables.  Let  the  error  in  the  predic¬ 
ted  variables  be  denoted  by  Z<M,y),  X(x,y) ,  and  Y(x,y) ,  where 
<x,y)  gives  the  spatial  location  and  it  is  assumed  that 
X(x,y)  »  k^Z^tx.y)  and  Y(x,y>  »  k^Z^tXjy).  The  subscripts  x  and 
y  denote  partial  dif f erentiation  with  respect  to  x  and  y 
respectively.  Assume  that  the  errors  in  the  predicted  values  are 
stationary  (that  is,  the  statistics  do  not  depend  on  (x(y)),  and 
have  zero  mean.  Using  EC.1  to  denote  the  expected  value,  or 
ensemble  average,  the  spatial  covariance  function  for  Z,  as  a 
function  of  "lags"  s  and  t,  is 

R (s,t)  «  ECZ (x ,y) Z (x+s,y+t) 1  -  ECZ (x-s,y-t> Z (x ,y> 1  . 

The  latter  equality  follows  from  stationari ty.  Under  the  assump¬ 
tion  that  the  order  of  partial  differentiation  and  the  expected 
value  can  be  interchanged,  the  cross  covariance  functions  and  the 
covariance  functions  for  the  derived  variables  are  found  in  the 
manner  illustrated  here.  Of  course,  it  is  assumed  throughout 
this  paper  that  the  necessary  derivatives  exist. 


ECZ(x,y)X(x+s,y+t)  3  ■  ItKN.ylk^Z^x+t.yH)]  ■ 
ECZ(x,y>fcZ<x-^*,y+t>3  -  kjECZ <x ,y > Z <x+*,y+t > 3#  -  k1Rs(*«t>  * 
Whll* 

ECX(x,y)Y<x+s,y+t>3  *  ECX <m  ,y> k2Zy  <x+s,y+t>  3  - 
ECX<K,y)k2Zt<N+«,y+t>3  -  kgECX  <x  ,y>  Z  <x+s,y+t>  3t  - 
kjECkjZ^ <x— s,y— t> Z  <x ,y>  3^  -  k^C-k^  <x-s,y-t >  Z  <x  ,y >  3 - 
-klk2ECZ<x,y)Z(x^«,y+t>  3t>  -  "klk2Rts <m,t>  . 

Not*  that  while  tha  covarianca  functions  ars  symmetric,  tha 
cross  covarianca  functions  ara  antisymmetric,  which  accounts  for 
tha  sign  change  that  coaas  froa  changing  tha  ordar  of  tha  product 
in  tha  expected  valua.  This  naans,  among  othar  things,  that  tha 
cross  covarianca  must  ba  zero  at  zero  lag  values.  This  behavior 
can  ba  seen  in  the  function  plats  in  Bsrgman  Cl 13  and  Schl attar, 
at.  al.  C 103 . 

2.2  Home  Necessary  Properties 

In  ordar  for  tha  covariances  of  tha  derived  functions  and 

tha  cross  covarianca  functions  to  exist,  certain  conditions  must 

ba  satisfied  by  the  function  R(s,t>.  Thase  have  been  alluded  to 

by  Buall  C173,  and  are  given  in  Julian  and  ThiSbaux  C183,  where 

R  <s>  R  <s> 

lie  -= -  is  finite,  and  limC  -= - R  (s)  3  -  0 

sX)  *  sX>  *  ** 

s  in  this  aquation  represents  tha  lag  distance  (lag  in  tha  above 
2  2  1/2 

was  (s  *  t  )  ),  and  R(s)  is  an  isotropic  covariance  function. 

Whan  ana  considers  that  R^tO)  must  be  zero,  the  first  limit  is 
tha  definition  of  tha  second  derivative  at  s"0,  hence  existence 
of  tha  limit  means  that  tha  covariance  function  must  be  twice 


dilfarantiabla  at  v<.  Tha  sac on d  limit  than  says  that  tha 
sacond  dar&vativa  is  continuous  at  s»0.  Thus  tha  thaoraa  given 
by  Julian  and  Thidbaux  can  ba  simplified! 

Thaoraa  li  If  R(s)  is  an  isotropic  covarianca  function  for  Z  in 
two  di mansions,  than  tha  covarianca  functions  for  tha 
partial  darivativas  of  Z  exist  at  s»0  if  and  only  if  R(s)  is 
twica  continuously  dif f arantiabla  at  s*0. 

2.3  Anisotropic  Functions 

It  has  baan  contandad  that  isotropic  covarianca  functions  do 
not  adaquataly  modal  tha  foracast  arror  statistics  and  that  gains 
can  ba  mada  by  using  anisotropic  functions.  Saa  Thidbaux 
C133 , C19D , C203 ,  and  Thidbaux,  at.  al .  C211  for  davalopmant  and 
discussion  of  product  forms  of  covarianca  functions.  Usa  of 
products  of  single  dimensional  functions  has  tha  advantage  of 
carrying  over  desirable  properties  to  higher  dimensions,  as  Mali 
as  being  able  to  use  essentially  one  dimensional  structures  and 
techniques.  On  the  other  hand,  perusal  of  contour  plots  of 
product  functions  show  that  zero  crossings  of  the  functions  occur 
along  grid  lines,  and  it  is  easy  to  see  this  will  always  happen. 
This  may  be  undesirable  behavior,  and  almost  certainly  it  is  not 
the  kind  of  behavior  seen  in  the  error  statistics. 

Another  form  of  anisotropy  is  possible,  one  which  results 
from  scaling  differently  in  two  orthogonal  directions,  then  using 
an  isotropic  functions  in  the  scaled  variables.  This  would 
result  in  the  zero  crossings  in  the  contour  plots  of  the  function 
being  ellipses  with  axes  in  the  two  directions,  and  all  contours 
having  the  same  shape.  The  eccentricity  of  the  ellipse  is  a 
measure  of  the  anisotropy  of  the  error  statistics.  It  would  be 
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May  to  allow  rotation  along  with  tho  scaling  to  obtain  ellipse 
of  constant  "distance"  with  any  axis  orisntation.  For  a  discus 


sion  of  this  type  of  anisotropic  correlations,  sss  Isaman  C223, 
and  Buoll  and  80 aaan  C233.  Ths  propsrtios  of  any  such  functions 
ara  thoss  of  isotropic  functions,  of  course,  sines  ths  anisotropy 
arises  purely  free  a  rotation  and  scaling. 

2.4  Isotropic  Functions 

The  use  of  isotropic  functions  in  two  or  sore  di sens ions 
that  have  been  derived  free  one  disensional  considerations  can 
possibly  lead  to  nonpositive  definite  functions.  For  example, 
Ripley  C243 (p  11)  quotes  a  result  of  Matern  C2S3,  which  gives  a 
lower  bound  for  isotropic  positive  definite  functions  in  several 
disensions.  The  result  ssans  that  positive  definite  functions  in 
two  dimensions  are  necessarily  bounded  below  by  -0.403  (the 

minimum  value  of  JQ(s)  ),  while  in  three  dimensions  the  bound  is 
-0.218  .  Thus  any  oscillatory  positive  definite  function  in  one 
dimension  that  takes  on  values  less  than  -0.403  cannot  be  an 
isotropic  positive  definite  function  in  two  dimensions.  A  posi¬ 
tive  definite  function  with  parameters  to  separately  control  the 
oscillation  frequency  and  the  decay  can  probably  be  made  into  a 
nonpositive  definite  isotropic  function  in  two  dimensions.  For 
example,  an  exponentially  damped  cosine  function,  f(s)  ■  cos (as) 
exp(-bs),  can  be  made  nonpositive  definite  by  suitable  choice  of 
parameters,  say  a  <■  S  and  b  ■  .1  .  This  result  also  applies  to 
other  candidates  for  Isotropic  correlation  function  models,  as 
will  be  shown  later. 

There  is  a  one-to-one  correspondence  between  covariance 


function*  In 


dia*n*lon  and  isotropic  covariance  functions  In 


too  dioan signs.  Using  tho  so  calisd  "turning  band"  sothod, 
Hathsron  C2A3  givss  a  way  of  ganaratlng  an  Isotropic  d-di son- 
si  anal  covarianca  function  fros  a  on#  disanslonal  covarianc# 


function.  Tho  relation  is 

rl 

Cd<s>  -  K  I  CjCvs)  <l-> 

J  A 


va,<d-3)/2  dv  f 


o  K  is  a  constant  that  is  uni sport ant  for  our  purposes.  In 


two  dieensions,  this  givss 


f1 

(S)  ■  K  Cx 


2  -1/2 

(vs) <l-v*>  4  *  dv  . 


It  is  possible  to  invert  Hathsron *s  relation  to  show  a  one-to-one 
relationship.  A  sketch  of  tha  inversion  process  follows. 
Esploying  a  change  of  variables  in  the  previous  expression  gives 


c2<s> 


-  K/c  i 


ctx.2-t2)-1/2  dt.. 


2  2 

then  eaking  further  change  of  variables,  s  «  x,  and  t  -  y,  yield 


C2<«‘/2> 


-  k/cCj 


(y1/2> <2y>”l/23 <y— x>“1/2  dt  . 


1 /2  —1  /? 

This  is  Abel's  equation  for  K  Cj(y  ><2y)  ,  and  the  solution 

is  well  known  (see  Hochstadt,  C273),  given  by 


V”*’  -  K  V'2  t  C2 


C,Cyl/2> <x-y>“l/2  dy  , 


tore  K'  is  a  different  constant.  Substitution  for  s  and  t  once 


again,  give* 


:i<«>  - K  - 

v  A 


2  2  -1/2 

<t>  <s  -t^>  4/^2t  dt  . 


The  correspondence  between  covariances  in  one  dimension  and 


th rmm  diwmiioni  is  sssisr  to  invert,  and  is  given  by  Ripley 
C243.  There  the  relation  is 


C^(s)  ■  I  C| (vs)  dv,  and  Cj <s>  ■  ^  CsC^(s)3  . 

J0 

While  this  characterization  of  eulti dimensional  isotropic 
covariance  functions  is  interesting,  and  can  in  fact  be  used  to 
generate  isotropic  eulti dimensional  covariance  functions,  it  does 
not  easily  answer  the  question  as  to  whether  or  not  a  particular 
one  dimensional  function  is  an  isotropic  positive  definite  func¬ 
tion  in  more  dimensions.  One  way  to  answer  such  a  question  is  to 
use  the  characterization  of  positive  definite  functions  as 
Fourier  transforms  of  probability  density  functions  (or  alterna¬ 
tively,  as  functions  whose  Fourier  transform  is  positive).  The 

Fourier  transform  of  an  isotropic  function  C(s>  in  two  dimensions 

1/2 

becomes  (essentially)  the  Hankel  transform  of  s  C(s).  It  may 
be  considerably  easier  to  look  at  the  one  dimensional  Fourier 
transform.  It  would  then  be  useful  to  have  a  sufficient  condi¬ 
tion  on  the  Fourier  transform  of  the  function  which  would 
guarantee  it  is  an  isotropic  positive  definite  function  in  two 
dimensions.  Such  a  condition  will  now  be  derived.  Let  Cjts)  be 
a  positive  definite  function  of  one  variable.  From  the  charact¬ 
erization  in  the  previous  section,  it  can  be  shown 


(2.1)  Cj(s) 


cos(rs) 


h(r)  dr. 


for  some  probability  density  function  h(r)  (i.e. ,  h(r)>0,  with 
integral  equal  1).  The  problem  is  then  to  determine  the  condi¬ 
tions  that  will  make  (s)  the  two  dimensional  Fourier  transform 
of  an  Isotropic  probability  density  function.  Such  a  transform 


is  nacMtarlly  isotropic.  A  -function  g(s)  is  sought  so  that 


(2.2) 


!i‘r»  “  / 

J  r% 


JA(rs)  s  g(s)  ds 


This  expression  is  invsrtsd  using  ths  Hsnksl  transform,  giving 


(2.3) 


g  (s) 


/. 


Jg(sr>  r  Cj (r>  dr. 


Than,  using  (2.1)  in  (2.3),  and  intarchanging  tha  ordar  of 
integration,  followed  by  integration  by  parts  yields 


g(s) 


Jq(s r)  r  (  /  cos(tr)  h(t>  dt>dr  ■ 


h(t)  <  /  Jq <ar )  r  cos(tr)  dr)  dt 


-  /  J0(sr)  r  (  J 

J0  Jo 

[It,  <  f 

J o  J  o 

/h<t.  ,-^f 

Jo 

«/; 

VA  A 


Jq(st)  sin(tr)  dr)  dt  *» 


<t)  <  )  Jq(st)  sin(tr)  dr  )  dt  , 


and  than, 

(2.4) 


g  (s) 


(t)/(t2-s2) 1/2  dt 


Tha  last  equality  uses  tha  Hankel  transform  of  r"*1/2sin <tr ) . 


Zn  order  for  g(s)  to  ba  a  probability  density  function  it 
must  be  nonnegative  with  integral  equal  to  one.  It  is  easy  to 
show  (again,  interchanging  the  order  of  integration)  that  the 
integral  is  equal  to  one.  It  is  more  difficult  to  show  necessary 
and  sufficient  conditions  for  the  nonnegativity  of  g(s).  The 
above  relations  summarized  givesi 


Theorem  2i  A  sufficient  condition  for  C.(s)  to  be  a  valid  iso¬ 
tropic  covariance  function  in  two  dimensions  is  that  h(t)  be 
a  monotone  decreasing  (h'(t)<0)  function. 
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This  condition 


unnscossoriiy  restrictive,  and  dlffl 


cult  to  uso  sinco  tho  condition  is  on  ths  Fourier  cosino 
transform  of  (s)  rather  than  C^  <s>  itself.  Nonetheless,  the 
condition  can  be  used  to  shou  the  following  interesting  results. 
1.  Consider  the  exponentially  daaped  cosine  function, 

C(s)  ■  coa  (ae)  e”***. 

The  Fourier  cosine  transfora  of  this  function  is 

bf.t  m  e (f )  ■  h H£±af!±&f!12  — 

h<t>  r  <C)  it)  Cba+(a-t>*3Cb*+<a+t>*3  . 

Inspection  of  h'(t)  shows  that  if  b2  £  3a2,  it  is  nonpositive  for 
all  t,  and  hence  h(t)  is  eonotone  decreasing  under  that 
constraint. 


II .  Consider  the  second  order  autoregressive  (SOAR)  covar¬ 


iance  function, 

C(s)  ■  Ceos (as)  +  (b/a)sin(as) 3e 


-bs 


The  Fourier  cosine  transform  of  this  function  is 


h(t)  ■  F(C>  <t)  - 


2b<b»+a»> 


Cb*+  <t-a) a3Cba+ C t+a) *3 


Inspection  of  h'(t)  reveals  that  if  b2  >  a2,  it  is  nonpositive 
for  all  t,  and  thus  h(t)  is  monotone  decreasing  under  that 
constraint.  We  see  that  each  of  the  above  C(s)  is  an  isotropic 
positive  definite  function,  hence  is  a  covariance  function  if  the 
appropriate  inequality  on  the  parameters  is  satisfied. 

ZZZ.  Consider  the  special  case  of  the  damped  cosine 


function 

C(s>  -  CA  +  (l-A) cos (as) 3/(1  +  (bs)2)172  . 
The  Fourier  transform  of  this  function  is 


h(t)  *  F(CXt)  - 

<2b>”1f  (1-A)CK0<  It-al/b)  K0<  lt+al /b>  3  +  AKQ(t/b>>  . 


Because  ttw  aodlflad  hual  function  KQ  boco— ■  unboundod  «•  tho 
arguMnt  tondo  to  xarof  for  A  f  1,  tho  Fourier  transfora  oust  bo 
increasing  os  t  approochos  •  through  vsluos  sosllor  than  a.  Far 
t  grootor  than  a,  and  possibly  for  soos  valuos  soallor  than  a, 
tho  function  is  docroasing.  Thus  tho  sufficient  conditions  givon 
abovo  aro  not  sot,  and  it  is  aasy  to  find  conf igurations  of  (x,y> 
points  and  paraostsr  valuos  A,  a,  and  b  for  Mhich  tho  resulting 
"covariance"  Matrix  is  not  positive  definite.  Tho  too  di son- 


si  onal  Fourior  transfors  of  C(s>  (tho  Hankel  transfora  of 
s1/2C(s>  )  has  thus  far  gone  unsolved,  so  it  is  prsssntly  unknown 


if  there  aro 


or  valuos  (ether  than  for  A  «  1)  that  will 


yield  a  positive  definite  function. 

XV.  Consider  the  Bessel  function  J^(as).  The  Four i mr 
transforo  of  this  function  is 


h  (t)  -  F  (C)  (t)  • 


This  function  is  easily 


0  ,  t<a 

t‘/2/<t2-.2>l/2  ,  t>.  . 


to  be  Monotone  decreasing  for  t>a. 


and  thus  the  Bessel  function  Jq(ss)  is  an  isotropic  covariance 
function  in  two  diaensions.  Application  of  this  relation 


requires  attention  to 


technical  details  because  of  the 


infinite  Jump  discontinuity  at  t«a. 

The  above  results  concerning  several  functions  proposed  for 
use  as  isotropic  covariance  functions  in  two  diaensions  is 
useful.  The  lack  of  results  and  enplrical  evidence  against  the 
daaped  cosine  being  positive  definite  negate  the  results  noted  in 
the  next  section  where  we  see  that  the  fitting  power  of  the 


the  next  section  where  we  see  that  the  fitting  power  of  the 
function  is  very  good.  These  aspects  of  the  function  will  be 


discussed  further  in  the  next  section. 


2.9 


y 


The  eontmti  of  this  Met&on  contain  aoao  u«*4ul  in-f  or  action 
for  tho  construction  of  isotropic  positive  definite  functions  and 
tasting  of  functions  for  positive  definiteness.  When  possible, 
the  two  diaensional  Fourier  transfora  of  Cj <s)  can  be  used  to 
decide  whether  or  not  the  function  is  positive  definite.  When 
the  two  diaensional  Fourier  transfora  cannot  be  obtained  in 
closed  fora.  Theorem  2  can  give  soae  inf or action  if  the  one 
diaensional  Fourier  transfora  is  available  in  closed  fora.  While 
the  sufficient  condition  given  by  Theorem  2  is  not  necessary,  it 
has  been  shown  to  be  useful  in  investigating  soae  functions  which 
have  been  proposed  for  use  as  isotropic  covariance  functions  in 
two  dimensions. 


3.0  8oae  Experiments  with  Isotropic  Covariance  Functions 


3. 1  Background 

The  work  reported  in  this  section  is  intended  to  help  deter¬ 
mine  something  about  the  overall  fitting  properties  of  various 
suggested  covariance  functions.  The  term  Hoverall  fitting  prop¬ 
erties'*  is  meant  to  include  not  only  the  ability  of  the  function 
to  model  a  reasonably  complicated  true  covariance  function,  but 
also  its  performance  when  used  in  a  statistical  interpolation 
scheme  with  several  different  observation  patterns. 


The  approach  for  this  project  was  to  begin  with  published 
data  from  an  actual  case,  and  then  construct  a  covariance  func¬ 


tion  using  a  least  squares  fit  to  the  data  from  a  certain  class 
of  covariance  functions.  This  model  is  used  to  define  the 


"truth"  model.  Functions  from  other  classes  were  then  fit  to  the 


*  *'•  ri  v 


mm  dtk«t  again  In  tha  Imt  squares  sense,  and  tha  performance 
of  thaaa  "aeeueed"  covariance  functions  ware  Matured  against 
that  of  tha  optiaua  aodel.  Tha  results  to  be  discussed  give  soas 
insight  into  what  classes  of  functions  have  adequate  fitting 
ability  for  aodel ing  actual  forecast  error  statistics,  and  also 
show  how  auch  skill  is  lost  (in  the  idealized  case)  by  use  of 
inaccurate  covariance  functions. 

The  results  given  here  consist  of  reprssentati ve  plots  of 
assueed  correlation  functions  together  with  the  correlation 
function  defined  as  "truth",  and  contour  plots  of  sqm  of  the 
resulting  expected  errors.  Tha  tables  show  expected  root— Man- 
square  (eras)  errors  (relative  to  the  standard  deviation  of  the 
background  error)  over  three  grids  of  paints  and  associated 
observation  locations.  The  expected  errors  were  coeputed  as  in 
Beaean  (1983).  The  results  obtained  with  various  assumed  corre¬ 
lation  functions  in  the  SI  scheM  are  discussed  in  detail. 
Additional  plots  are  given  in  Franks  C283. 

3.2  The  Model  Correlation  Function 

The  data  for  the  covariance  function  was  obtained  (by  hand) 

from  Lonnberg  C293.  The  data  taken  was  plotted  points  from  a 

covariance  function  of  the  type  used  by  European  Canter  for 

Medium-Range  Weather  Forecasts,  in  this  case  a  five  term  (l.e., 

n*S)  Bessel  series  of  the  fore 
n 

(3.1)  y^Ai<TQ(s»kt/R)  ♦  AQ  , 
i-1 

where  k^  is  the  ith  zero  of  the  Bessel  function  Jq(s),  and  R  is 
the  radius  of  the  region  of  interest.  This  function  is  positive 


tfaflnlta  as  an  Isotropic  function  in  two  dimensions  provided  tho 


coefficient*  ere  ell  positive.  In  Lonnberg,  R  wee  2000  k«. 

In  this  work,  distance  was  eeasured  in  degrees,  and  the  radius 
was  scaled  to  30°. 

Least  square*  fits  to  the  data  by  functions  of  the  type 
(3.1)  for  four,  five,  and  six-terms  were  computed.  While  the 
original  paper  indicated  a  series  with  five  terms  generated  the 
data,  it  was  found  that  six-terms  yielded  all  positive  coeffi¬ 
cients  and  a  significant  reduction  in  the  residual  over  five 
terms.  Thus,  it  was  decided  to  adopt  the  six  term  series  as  the 
"truth"  covariance  function.  This  six  term  series  would  also  be 
marginally  harder  to  approximate  using  other  classes  of 
covariance  functions.  The  data  and  the  fits  using  four  and  six 
terms  are  shown  in  Figure  1,  and  the  coefficients  are  given  in 
Table  1,  along  with  other  data.  The  intercept  values  of  the 
approximations  were  0.8270  and  0.8592  for  four  and  six  terms, 
respectively.  This  occurs  because  the  data  represents  the 
spatial  correlation  of  the  background  plus  observation  error, 
thus  the  intercept  is  a  function  of  the  ratio  of  the  standard 
deviations  of  background  and  observation  error.  The  effects  of 
this  kind  of  discrepancy  will  be  discussed  in  section  3.3.  The 
correlation  function  for  background  error  is  the  approximation 
normalized  to  have  value  one  at  s«0,  of  course. 

3.3  The  Brid  and  Observation  Point  Sets 

Three  grids  and  associated  point  sets  were  selected  for 
studying  the  expected  errors  of  statistical  interpolation  schemes 
based  on  various  assumed  covariance  functions.  All  were  based  on 
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til*  approMlMta  location*  of  radiosonde  data  (from  Wahba  and 
Wendelberger  C303  and  Shil,  ot.al.  C313)  within  the  aoloctod 
grid.  Each  grid  covorad  a  region  that  was  30°  in  longitude  and 
20°  in  latitude,  and  the  three  Here  chosen  to  represent  a  dense 
observation  set,  a  partially  dense  observation  set,  and  a  sparse 
observation  set.  The  regions  correspond  to  the  middle  United 
States  nith  36  observations,  the  eastern  United  States  and  west¬ 
ern  Atlantic  Ocean  with  25  observations,  and  the  middle  Atlantic 
Ocean  with  3  observations.  For  reference  purposes,  the  three 
regions  will  be  referred  to  as  the  MU8  (Mid-US),  EC  (East  Coast), 
and  MA  (Mid-Atlantic)  regions.  The  regions  and  the  observation 
locations  can  be  seen  in  Figures  2-3,  parts  (b) ,  (c) ,  and  (d) , 
respectively.  The  regions  were  gridded  at  2.5°  intervals  for 
purposes  of  computing  expected  errors,  although  the  erms  errors 
given  in  Table  1  are  only  over  the  interior  grid  points  to  mini¬ 
mize  edge  effects.  Use  of  interior  grid  points  for  this  purpose 
is  valid  since  on  a  sphere  it  is  not  necessary  to  interpolate  to 
the  boundary  points.  For  contouring  purposes  the  fields  were 
interpolated  to  finer  grids  using  bicubic  spline  interpolation. 

3.4  The  Assueed  Correlation  Functions  and  Results 

The  families  of  assumed  correlation  functions  fell  into  5 
classesi  (1)  Bessel  function,  (ii)  negative  squared 
exponential  (sometimes  called  Gaussian),  (iii)  autoregressi ve  of 
order  two,  (iv)  autoregressi ve  of  order  three,  and  (v)  damped 
cosine.  They  will  be  discussed  in  turn,  along  with  the  results. 
Plots  of  the  assumed  correlation  functions,  along  with  the 
"truth"  correlation  function,  are  shown  in  part  (a)  of  Figs.  2-3. 


Far  fitting  purposes*  each  included  a  multiplicative  paraaatar 
that  determined  tha  i«0  intercept,  and  was  aubaaquantly  droppad 
to  obtain  tha  corralation  function.  Tha  valua  of  thia  paraaatar 
ia  of  intarest,  however ,  bacauaa  dropping  it  ahifta  tha  curva 
(upward)  to  paaa  through  tha  point  (0,1),  and  thua  diffarant  fita 
may  ba  ahiftad  by  diffaring  amount*,  which  ultimataly  affacta  tha 
fit  to  tha  background  arror  corralation  function. 

Racall  that  tha  arma  arrora  givan  in  Tabla  1  ara  given  aa  a 
fraction  of  atandard  deviation  of  tha  background  arror.  Tha 
ratio  of  tha  atandard  daviation  of  tha  observation  errors  to  the 
standard  deviation  of  tha  background  arrora  was  1/3. 

(i)  Bessel  Function 

Tha  reference  expected  errors  ware  computed  using  the  actual 
correlation  function  model,  given  by  Eq.  (3.1)  with  coefficients 
as  givan  in  Table  1.  The  results  are  given  in  Table  1,  and  are 
tha  smallest  expected  errors  that  can  be  obtained  using  a  correc¬ 
tion  to  background  scheme,  that  is,  they  are  truly  optimum.  The 
correlation  function  is  shown  in  part  (a)  of  Fig.  2,  while  the 
contour  plots  of  the  expected  error  for  each  of  the  three  grid/ 
observation  sets  is  shown  in  parts  (b),  (c) ,  and  (d) . 

Tha  results  using  a  four  term  Bassal  function  ara  givan  in 
Tabla  1.  Because  tha  intarcapt  of  tha  fit  to  tha  data  is  0.8270 
versus  0.8592,  normalization  to  valua  one  produces  a  curva  that 
is  predominately  above  that  for  tha  modal  correlation  function, 
especially  for  small  distances.  Tha  result  of  the  poor  approxi¬ 
mation  for  small  distances  is  most  pronounced  over  MU8.  The 
effect  was  small  over  the  sparse  part  of  EC  and  over  MA. 


(i i )  Negative  Squared  Exponential  (NSE) 


The  NSE  has  bssn  recognized  aa  i nad squats  for  sodsling  arror 
covariancss  for  boss  time,  and  tha  rsaults  obtalnsd  Kara  confirs 
that.  Tha  assusad  fors  of  tha  function  was 

(3.2)  A  +  <l-A)a"<B/b>a  . 

This  function  is  positive  dafinita  as  an  isotropic  function  in 
tso  diaensions  for  0< A< 1 . 

Tha  initial  fit  Nas  not  obtained  by  least  squares,  but 
siaply  by  attempting  to  fit  tha  modal  correlation  reasonably  wall 
for  small  distances,  taking  A-0.  Tha  fit  is  reasonably  good  up 
to  about  6°,  and  quits  poor  at  graatar  distances.  Tha  expected 
errors  are  similar  in  magnitude  to  tha  expected  errors  for  tha 
four  tars  Bessel  function,  axcapt  over  MA,  where  the  errors  are 
larger.  However,  since  the  errors  over  MA  tend  to  be  large 
anyway,  tha  relative  effect  is  not  as  great  as  one  mig*>t  expect. 

The  second  attempt  was  by  least  squares  for  the  parameters  A 

and  b.  Because  the  NSE  is  too  flat  near  the  origin,  this  process 

yielded  an  intercept  value  of  0.8060,  shifting  the  correlation 
function  so  that  it  is  entirely  above  the  model  correlation 
curve.  This  results  in  even  poorer  perf ormanca  over  MUS  and  EC 
than  the  previous  modsl ,  dus  to  ths  inaccurate  representation  for 
small  distances.  The  performance  over  HA  was  bettsr  than  the 
above. 

Due  to  the  poor  performance  (compared  to  the  above)  obtained 
by  adding  a  constant  to  the  basic  NSE  it  was  decided  to  attempt 

to  find  a  better  fit  by  trial  and  error.  No  claim  is  made  about 

any  optimality  for  this  function.  The  results  in  Table  1  demon¬ 
strate  that  it  is  probably  not  possible  to  obtain  good  results 


overall  with  a  function  from  tho  N8C  family,  and  certainly  not 
for  tho  promo nt  oodol  corrolation  function. 

(iii)  Autoregressive,  Order  Two  (SOAR) 

Tho  SOAR  oodol  has  boon  suggested  as  appropriate  by  Yudin 
C323,  ThiObaux  C 133  and  this  is  supported  by  simulations  due  to 
Balgovind,  et.al.  C333.  This  is  the  model  that  is  being  incor — 
p orated  into  the  U.  8.  Navy  MNP  models.  The  formula  given  here 
includes  a  constant  term  which  is  not  part  of  the  SOAR  model,  but 
which  has  been  noted  to  improve  performance  considerably 
(ThiObaux,  et.al.,  C211),  and  those  results  are  confirmed  here. 
The  SOAR  function  with  additive  constant  is 

(3.3)  A  +  (1-A) Ceos (as)  +  (b/a)sin (as) 3e  b* 

This  function  is  positive  definite  (in  two  dimensions)  whenever 
a<b,  and  0< A< 1 .  In  all  cases  investigated  here,  and  as  has  been 

reported  elsewhere,  (e.g.,  ThiObaux,  et.al.,  C213),  the  parameter 
a  tends  to  be  essentially  zero.  In  this  case  the  function  re¬ 
duces  to 

(3.3a)  A  +  (1-A) Cl  ♦  bs3e~b*  . 

The  initial  attempt  was  a  least  squares  fit  to  the  data  with 
A«0.  The  intercept  obtained  was  0.7977,  with  the  resulting 
correlation  curve  then  being  considerably  above  the  model  corre¬ 
lation  curve  between  0°  and  15°.  The  performance  was  only 
slightly  better  than  with  any  of  the  previous  correlation 
functions.  It  was  then  decided  to  attempt  a  least  squares  fit 
with  the  intercept  constrained  to  be  0.8S92,  the  same  as  obtained 
for  the  model  correlation  function,  but  again  with  A»0.  Table  1 
shows  marginal  improvement  for  all  three  grid/observation 


p *t tarns.  A  third  attempt  included  A  in  ths  lssst  squarss  fit, 

with  no  constraint.  This  rasultad  in  a  such  closer  eatch  to  the 
aodel  correlation  function,  although  the  Intercept  of  0.6441 
eoved  the  assumed  correlation  curve  above  the  model  curve  for 
much  of  the  interval.  The  fit  and  resulting  expected  error 
contours  are  shown  in  Figure  3.  Table  1  shows  thsre  is  consid¬ 
erable  improvement  over  all  previous  results,  the  most  improve¬ 
ment  being  for  MUS,  and  the  least  for  MA. 

(iv)  Autoregressive,  Third  Order  (TQAR) 

The  use  of  the  TOAR  model  has  been  investigated  by  ThiSbaux, 
et  al.  C213,  including  an  additive  constant.  The  formula  is 

(3.4)  A  +  (1-A) C (acos(as)  +  bsin(as))e  b*  +■  ce  c*3  , 

where  the  coefficients  a,  b,  and  c  are  functions  of  a,  b,  and  c 
given  by 

a  -  <3b2-a2-c2)ac/D  ,  b  -  (b2-3a2-c2)bc/D  , 
c  -  -2  <b2+a2>  ab/D  ,  where  D  -  (3b2-a2-c2)ac-2(b2+a2)ab  . 

It  is  unknown  what  restrictions  (beyond  0<A<1)  on  the  parameters 
are  required  to  ensure  the  function  is  positive  definite  as  an 
isotropic  function  in  two  dimensions. 

The  data  was  fit  by  least  squares  with  the  TQAR  function 
(3.4).  The  intercept  was  0.8651,  which  resulted  in  the  curve 
being  slightly  below  the  model  correlation  curve  over  most  of  the 
range.  Overall,  the  fit  was  quite  close  and  better  than  any  of 
the  previously  discussed  functions.  The  results  in  Table  1  show 
very  close  agreement  with  the  optimum  possible  for  all  three  of 
the  grid/observation  sets. 


(v)  Diapid  Cosin* 


Ths  damped  cosins  function  hss  boon  suggested  by  Thlbbaux 
C193  and  Seaman  and  Hutchinson  C343.  Tho  for aula  is 
(3.3)  CA  +  (l-A)eos(as) 3/Cl  +  (bs)23c  . 

It  is  unknown  whsthsr  tho  function  is  positivo  dofinlts  as  an 
isotropic  function  in  two  dimensions,  but  tho  svldonco  in  soction 
2.4  (for  c*0.3) ,  while  inconclusive,  seems  to  lndlcata  it  Is  not. 
In  practice,  of  course,  the  function  may  bo  positivo  definite 
whan  tho  observation  points  ara  restricted  to  certain  regions. 

The  data  was  fit  with  the  function  (3.5),  under  the  restriction 
c^O.S.  The  intercept  was  0.8363,  which  resulted  in  a  very  slight 
raising  of  the  curve  relative  to  the  model  correlation  function. 
Tha  resulting  fit  is  excellent  for  small  distances  and  very  good 
over  the  entire  range.  Table  1  shows  that  this  function  gives 
tha  best  results  of  all  the  functions  testad. 

(vi)  Variations 

The  expected  error  computations  for  a  number  of  variations 
of  the  above  functions  were  also  performed.  Tha  principal 
variation  was  to  fit  the  data  only  over  the  first  half  of  the 
interval,  (0°,13°).  The  effect  of  this  was  to  generally  (though 
not  always)  increasa  the  ores  errors  over  MU8  and  EC,  while  not 
effecting  the  results  over  HA.  In  the  damped  cosine,  the 
exponent  c  was  chosen  by  least  squares,  along  with  the  other 
parameters,  and  resulted  in  a  slightly  better  fit  to  the  correla¬ 
tion  function,  especially  at  larger  distances.  However,  the 
coefficient  A  was  slightly  greater  than  1.  Whatever  the  positive 
definiteness  properties  of  the  function,  having  A>1  will  certain¬ 
ly  make  it  non-positive  definite.  Although  no  graphical  results 


arm  shown,  the  coefficients  and  mrmm  srrors  ars  given  in  Table  1 
for  tha  additional  asauaod  corralation  functions. 

3.5  Mansi tl vi ty  of  tha  80AR  Modal  to  Paraaatar  Mi sspacifi cation 

In  ardor  to  dataraina  aora  coaplataly  tha  charactaristics  of 
tha  SOAR  modal,  soma  additional  calculations  wars  mada  to  deter¬ 
mine  tha  affact  of  ai sspacifi cation  of  tha  paramatars  in  tha 
corralation  function  or  tha  ratio  of  tha  standard  daviations  of 
tha  obsarvad  and  background  arror.  Tha  rasults  can  ba  mummed  up 
rathar  quicklyi  Tha  schaaa  is  mostly  insansitiva  to  such 
variations.  Figura  4  shows  a  family  of  4  corralation  functions, 
#4  baing  tha  SOAR  plus  constant  discussad  in  tha  saction  3.4, 
with  tha  othars  having  smallar  corralations  at  a  givan  distance. 
Figura  S  shows  tha  expected  RMS  arror  for  aach  of  tha  four  as  tha 
"assumed**  corralation,  whan  tha  "trua"  corralation  function  is 
•4.  With  tha  exemption  of  tha  sparse  MA  grid,  tha  expected 
errors  are  relatively  stable  under  significant  perturbations. 
Figura  6  shows  tha  sensitivity  to  tha  assumed  arror  ratio,  and 
ones  again,  it  is  observed  that  tha  expected  errors  are  quite 
stable. 

4.0  Conclusions 

Tha  principal  conclusion  to  ba  drawn  is  that  tha  correlation 
family  used  in  practical  analysis  should  embody  a  sufficient 
number  of  parameters  to  fit  tha  forecast  error  statistics  reason¬ 
ably  wall.  Further,  it  is  most  important  that  tha  data  be  fit 
accurately  for  small  distances.  In  order  to  ensure  a  batter  fit 
for  small  distances,  it  may  ba  worthwhile  to  enforce  tha  inter- 
capt  of  tha  corralation  for  tha  background  plus  observation 


error*  JJL  the  ratio  of  standard  deviations  of  the  two  errors  is 
known  y-  The  effect  of  scaling  to  obtain  the  correla¬ 

tion  function*  and  the  mppmrmnt  shift  up  or  down  can  possibly  be 
coapensated  for  by  artificially  varying  the  ratio  of  background 
to  observation  errors,  as  well,  although  it  seems  more  desirable 
to  enforce  this  ratio  in  the  correlation  function  fitting 
process. 

As  noted  above,  clearly  the  most  important  region  for  the 
fit  to  the  correlation  function  to  be  accurate  is  for  small 
distances.  Over  the  sparsely  observed  region,  MA,  and  to  a 
lesser  extent  over  the  EC  region,  the  overall  arms  errors  were 
only  slightly  affected  by  the  assumed  correlation  function.  In 
the  case  of  the  HA  region  it  is  noted  that  the  error  contours  are 
relatively  unaffected  except  near  the  observations.  Since  the 
errors  in  the  remote  part  of  the  region  dominate  the  overall 
error,  the  choice  of  assumed  correlation  function  has  relatively 
small  influence.  On  the  other  hand,  over  the  densely  observed 
region,  an  accurate  fit  at  small  distances  was  most  important. 

The  NSE  correlation  function,  while  not  performing  well,  illus¬ 
trates  the  above  nicely.  For  the  first  NSE  entry  in  Table  1, 
even  though  the  fit  is  poor  for  distances  of  more  than  6°,  the 
erms  errors  over  MUS  and  EC  are  smaller  than  the  best  fit  (last 
N8E  entry)  due  to  the  more  accurate  fit  for  small  distances  by 
the  former  function.  Of  course  the  erms  errors  over  MA  are 
poorer  for  the  first  case  due  to  the  very  bad  fit  at  large 
distances. 

There  appear  to  be  several  good  candidates  for  use  as  two 
dimensional  isotropic  correlation  functions,  including  SOAR, 


TOAR,  and  dupwi  cosine,  given  by  Eqs.  <3.3),  (3.4),  and  (3.3), 
raapactivaly.  While  tha  fitting  power  for  tha  1 attar  two  ara 
graatar  (thara  ara  a  graatar  nuabar  of  paraaatara  for  thoaa  two) , 
tha  choica  of  SOAR  aaaaa  raasonabla  and  adaquata  for  a  nuabar  of 
raaaonai  (1)  The  SOAR  (with  tha  additive  constant)  aabodlas  a 
sufficient  number  of  parameters  to  allow  oscillation  and  decay 
with  distance.  (2)  Tha  SOAR  has  soma  credibility  as  tha  spatial 
correlation  function  of  an  innovation  process.  However,  results 
are  for  one  dimension  rather  than  two,  an cap t  for  tha  results 
cited  previously  in  Balgovind,  at.al.  C333.  (3)  Tha  SOAR  was 

demonstrated  hera  to  be  positive  definite  as  an  isotropic 
function  in  two  dimensions,  undar  a  mild  restriction  on  tha 
parameters.  <4)  Whila  tha  TOAR  is  also  the  spatial  correlation 
function  (again  in  one  dimansion)  for  an  innovation  process, 
based  on  this  limited  study  it  does  not  appear  to  be  signifi¬ 
cantly  batter  than  SOAR.  (5)  Tha  positive  definiteness 
properties  of  the  TOAR  are  not  known,  although  it  is  certainly 
positive  definite  as  an  isotropic  function  in  two  dimensions 
under  some  restrictions  on  the  parameters.  (4)  Although  the 
fitting  ability  of  the  damped  cosine  seems  to  be  at  least  as  good 
as  the  TOAR,  and  it  is  positive  definite  in  one  dimension,  evi¬ 
dence  indicates  it  may  not  be  positive  definite  as  an  isotropic 
function  in  two  dimensions,  regardless  of  parameter  restrictions. 
The  availability  of  other  acceptable  alternatives  seems  to  make 
it  prudent  to  preclude  the  use  of  the  damped  cosine  in  practical 
situations. 

It  is  painted  out  that  all  of  the  functions  except  the  four 
term  Bessel  function  and  the  NBE  perform  very  well.  Table  i 


shows,  for  maiipla,  that  th*  SOAR  is  only  a  llttlo  mo re  than  IX 
of  tho  standard  daviation  of  tha  background  arror  poorer  than 
optlaal  over  MU8  and  EC,  and  1 aaa  than  .IX  poorer  over  MA. 

Finally,  it  is  noted  that  Mithin  the  SOAR  family,  81  ia 
quite  insensitive  to  eisspecif  i  cat  ion  of  the  correlation  paramo- 
ters,  even  to  an  extent  such  that  the  correspondence  would  appear 
to  be  much  less  between  two  members  of  the  family  than  between  it 
and  a  fit  by  the  NBE.  Thus  it  could  as  important  to  choose  the 
correct  family  of  correlation  functions  as  well  as  to  model 
properly  within  that  family.  In  addition,  ml sspecifi cation  of 
the  ratio  of  standard  deviations  of  the  background  and  observa¬ 
tion  errors  has  a  rather  small  effect  on  the  skill  of  the  method. 

This  work  has  focused  only  on  the  univariate  problem, 
whereas  in  practice  such  schemes  are  applied  to  the  multivariate 
one.  Further  work  is  necessary  to  determine  whether  the  nice 
results  obtained  here  carry  over  to  the  multivariate  case.  A 
further  investigation  of  the  effect  of  wind  observations  on  the 
analysis  of  pressure  height  and  wind  fields  is  anticipated. 
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PiTMity  Valu ••  and  araa  Errors 
Modal  Correlation  ia  81m  Tara  Bessel 


AmumnJ 


Correlation 

Function 

Paraeeters 

afb,c  A,Ai 

HUS 

ere s 

EC 

MA 

6  T  Bessel 

0.2474 

0.3333 

0. 1844 

0. 1031 
0.0362 
0.0554 
0.0400 

0.2667 

0. 3752 

0.7483 

4  T  Bessel 

0.2811 

0.3090 

0.2213 

0.0930 

0.0936 

0.3046 

0.4088 

0.7303 

NSE 

10.0 

0.0 

0. 3047 

0.4184 

0.7822 

N8E 

14.88 

0.3200 

0.3688 

0.3282 

0.7631 

NSE 

10.0 

0.2500 

0.3098 

0.4158 

0.7341 

SOAR 

0.0 

0. 1213 

0.0 

0.3034 

0.4022 

0.7562 

SOAR 

0.0 

0. 1374 

0.0 

0.2931 

0.3968 

0.7562 

SOAR 

0.0 

0.2033 

0.2722 

0. 2780 

0. 3839 

0.7491 

TOAR 

0.4732 
0 . 3828 
0.0914 

0. 1974 

0.2717 

0.3794 

0.7483 

Depd  Cos 

0.4749 
0. 1367 
O. 3000 

0. 9392 

0.2686 

0.3779 

0.7486 

N8E 

13.0 

0.0 

0.3619 

0.4414 

0.7649 

NSE 

12.31 

0.3205 

0.3474 

0. 4299 

0. 7393 

SOAR 

ja 

0.0 

0.2634 

0.3758 

0.2743 

0.3825 

0. 7493 

TOAR 

ja 

0.4468 
0. 1482 
0.0032 

-5.9965 

0. 2697 

0.3801 

0.7514 

Depd  Cos 

1.2236 
0. 1307 
0.3000 

1 . 0027 

0.2749 

0.3734 

0.7491 

Depd  Cos 

ja 

0.7009 
0.2069 
0. 3753 

1.0105 

0.2692 

0. 3779 

0. 7484 

Depd  Cos 

0. 79S7 
0.2350 
0.3317 

1.0147 

0.2706 

0.3784 

0. 7485 

Table  1 

*  These  correlation  functions  were  obtained  by  least  squares 
fit  over  the  interval  <0°,13°> 


FIGURE  CAPTIONS 


Figure  It  The  data  points  and  laaat  squares  fits  by  four  and  six 
tars  Bassal  functions. 

Figure  2i  (a)  Six  tarn  Bassal  corralation  function  (trua  and 
assumad).  (b)  Expactad  root-mean— square  arror  contours  for  tha  MUS 
grid  and  obsarvation  point  sat  for  tha  corralation  functions  shown 
in  (a).  (c)  As  in  (b)  for  tha  EC  grid  and  observation  point  sat. 

(d)  As  in  (b)  for  tha  MA  grid  and  obsarvation  point  sat. 

Figure  3»  (a)  Six  term  Bassal  correlation  function  (trua)  and 

second  order  autoregressive  plus  constant  correlation  function 
(assumad).  (b> ,  (c)  ,  and  (d)  As  in  corresponding  parts  of  Fig.  2 
for  tha  corralation  functions  shown  in  (a). 

Figure  4*  Four  2nd  order  autoregrassi ve  corralation  functions, 
as  in  Eq.  3.3a,  with  <b,A>  valuasi  #1  -  (0. 9,0.0))  #2  - 

<0.3, 0.0))  #3  -  (0.3-0.15))  #4  -  (0.2055,0.2722). 

Figure  5t  Expactad  RMS  errors  whan  tha  "trua"  correlation  is 
function  #4  with  various  assumad  correlation  functions  for  each 
of  tha  three  grid/observatian  sets. 

Figure  6t  Expected  RMS  errors  whan  the  assumed  ratio  of  the 
standard  deviations  of  the  observed  to  background  arror  is 
varied.  Actual  arror  ratio  is  1/3. 


LONNBERG  POINTS,  BESSEL  FTNS 


•  T» 


DISTANCE  (DEGREES) 


(•) 


EXPECTED  ERROR 

TRU=SIX  TERM  BESSEL  OBS/P3  =  0.333 
ASU=SX  TERM  BESSEL  OBS/TC  =  0.333 
9:33:04  4/2 V& 


(b) 

FIGURE  2 

34 


■n—fnwtnnri  rr 


EXPECTED  ERROR 

TRU-SIX  TERM  BESSEL  OBS/FG  -  0.333 

RSU-SIX  TERM  BESSEL  OBS/FG  -  0.333 

9:37:59  4/21/07 


<c> 


EXPECTED  ERROR 

TRU-SIX  TERM  BESSEL  OBS/FG  -  0.333 
flSU-SIX  TERM  BESSEL  OBS/FG  -  0.333 


>—■  ».  -  |  ■  ■  |  —  '  |—  ■■  ■  '  -  ■  f 

-«.o  K.t  -io.o  -33.0  -30.0  -a.a  -io. 

LONGITUDE 


cd> 

FIGURE  2 


TRUE  &  ASSUMED  CORR  FTNS 


s 


4 - 1  —  -  '  I - 1 - 1 - 1 - 1 

«  •  1C  19  10  >3  M 

DISTANCE  (DEGREES) 


<•> 


EXPECTED  ERROR 

TRU=StX  TERM  BESSEL  0BS/P3  =  0.333 
AS=£0  0.206  0.272  OBS/fc  =  0.333 
9:55:26  4/21/87 


-11X4 


—1074 


-174 

LONGITUDE 


-*24 


-124 


<b> 

FIGURE  3 


3A 


AR2  CORRELATION  FUNCTIONS 


DISTANCE  (DEGREES) 


FIGURE  4 


SENSITIVITY  OF  SI  TO  CORR  FTN 


-4 


’O 


» "  — —  I  1  T 

i  »  a 

ASSUMED  CORRELATION  FUNCTION  # 


T 

4 


FIGURE  S 


SENSITIVITY  OF  SI  TO  ERROR  RATIO 


FIGURE  6 


3* 


INITIAL  DISTRIBUTION  LIST 


A8ST.  FOR  ENV.  SCIENCES 
ASST.  SEC.  OF  THE  NAVY  (RStD) 
ROOM  3E731,  THE  PENTAGON 
WASHINGTON,  DC  20350 


OFFICE  OF  NAVAL  RESEARCH 
CODE  422AT 

ARLINGTON,  VA  22217 

COMMANDING  OFFICER 
FLENUMOCEANCEN 
MONTEREY,  CA  93943 

DR.  EDWARD  BARKER 
NAVAL  ENVIRONMENTAL  PREDICTION 
RESEARCH  FACILITY 
MONTEREY,  CA  93943 

DIRECTOR  (2) 

DEFENSE  TECH.  INFORMATION 
CENTER,  CAMERON  STATION 
ALEXANDRIA,  VA  22314 

Director  of  Research  Admin. 

CODE  012 

NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


CHIEF  OF  NAVAL  RESEARCH  (2) 
LIBRARY  SERVICE8,  CODE  784 
BALL8TQN  TOWER  «1 
800  QUINCY  ST. 

ARLINGTON,  VA  22217 

LIBRARY,  FLEET  NUMERICAL  (2) 
OCEANOGRAPHY  CENTER 
MONTEREY,  CA  93943 

LIBRARY  (2) 

NAVAL  POSTGRADUATE  8CHOOL 
MONTEREY,  CA  93943 

DR.  JAMES  GQER88 
NAVAL  ENVIRONMENTAL  PREDICTION 
RESEARCH  FACILITY 
MONTEREY,  CA  93943 

SUPERINTENDENT 
LIBRARY  REPORTS 
U.  S.  NAVAL  ACADEMY 
ANNAPOLIS,  MD  21402 

DEPT.  OF  MATHEMATICS 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


PROFESSOR  H.  FREDR I CK8EN 
CHAIRMAN,  DEPT  OF  MATHEMATICS 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 

DR.  RICHARD  LAU 
OFFICE  OF  NAVAL  RESEARCH 
800  QUINCY  8T. 

ARLINGTON,  VA  22217 


PR0FE880R  R.  FRANKE,  53F«  (10) 
DEPT.  OF  MATHEMATICS 
NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 

DR.  PAUL  F.  TWITCHELL 
OFFICE  OF  NAVAL  RESEARCH 
800  QUINCY  ST. 

ARLINGTON,  VA  22217 


DR.  THOMAS  RQ8M0ND 
NAVAL  ENVIRONMENTAL  PREDICTION 
RESEARCH  FACILITY 
MONTEREY,  CA  93943 

DR.  R.S.  SEAMAN 
AUSTRALIAN  NUMERICAL 
METEOROLOGY  RESEARCH  CENTRE 
P.O.  BOX  30S9AA 
MELBOURNE,  VICTORIA, 

AUSTRALIA,  3001 


DR.  H.  JEAN  THIEBAUX 

NATIONAL  METEOROLOGICAL  CENTER 

W/NMC2 

WWB 

WASHINGTON,  DC  20233 

Center  for  Naval  Analyses 
2000  N.  Beauregard  St. 

Alexandria,  VA  22311 


